function [] = diagnostics(pars,M_full,W_full,Nu,verbose)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vector of fixed parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rho                                     = NaN(numel(Nu.fixed_pars),1);
for p=1:numel(Nu.fixed_pars)
        Rho(p)                                  = pars.(Nu.fixed_pars{p});
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vector of targeted parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Xi                                      = NaN(numel(Nu.targeted_pars),1);
for p=1:numel(Nu.targeted_pars)
        Xi(p)                                  = pars.(Nu.targeted_pars{p});
end 
            
%%%%%%%%%%%%%
% Diagnostics
%%%%%%%%%%%%%

%Nmoments                                = size(M_full,1);
Nmoments                                = numel(Nu.targeted_moments);
Npars                                   = numel(Nu.targeted_pars);

% Cov matrix
dGdpars                                 = NaN(Nmoments,Npars);
h                                       = Nu.JacobianStepSize;
parfor n=1:Nmoments
    fun                                     = @(x) G_fun( x,Rho,Nu,M_full,n );
    df                                      = num_grad(fun, Xi, h);
    dGdpars(n,:)                            = df';
end    
Omega                                   = ( 1 + 1/Nu.Nsim )*( ( ( dGdpars' )*( W_full(Nu.log_targeted_moments,Nu.log_targeted_moments) )*( dGdpars ) )^(-1) ) ;
sd                                      = 1/sqrt(Npars)*sqrt(diag(Omega));
%sd                                      = 0.01*ones(Npars,1);

% OID test statistic
ObjTargMoms                             = SMM_obj( Xi, Rho, Nu, M_full, W_full, Nu.log_targeted_moments, 0 );
ObjAllMoms                              = SMM_obj( Xi, Rho, Nu, M_full, W_full, true(size(Nu.log_targeted_moments)), 0 );
if Nmoments > Npars
    OID                                     = 1/sqrt(Npars)*Nu.Nsim/( 1 + Nu.Nsim )*ObjAllMoms;
    dfr                                     = Nmoments-Npars;
    pval                                    = 1-chi2cdf(OID,dfr);
else
    [OID,dfr,pval]                          = deal(NaN);
end

% Monte-Carlo test for whether the estimated and empirical moments are different
Nmoments                                = size(M_full,1);
Nu.Nreps                                = Nu.Nreps_diags;
Nu.Seeds                                = Nu.Seeds_diags;
Nu.eps                                  = Nu.eps_diags;
OutputMoments                           = produce_moments(Xi,Rho,Nu);
MomDraws                                = NaN(Nu.Nreps,Nmoments);
MomStats                                = NaN(Nmoments,6);
for n=1:Nu.Nreps
    MomDraws(n,:)                           = (OutputMoments.M(:,n))';
end
for n=1:Nmoments
    MomStats(n,1)                           = M_full(n);
    MomStats(n,2)                           = mean(MomDraws(:,n));
    MomStats(n,3)                           = std(MomDraws(:,n));
    if MomStats(n,1) >= MomStats(n,2)
        MomStats(n,4)                           = sum( MomDraws(:,n) > MomStats(n,1) )/Nu.Nreps;
    else
        MomStats(n,4)                           = sum( MomDraws(:,n) < MomStats(n,1) )/Nu.Nreps;
    end
    MomStats(n,5)                           = prctile(MomDraws(:,n),5);
    MomStats(n,6)                           = prctile(MomDraws(:,n),95);
end

%%%%%%%%%
% Output
%%%%%%%%%    

if verbose 

    % Table 6
    Table6                                  = table(num2str([Xi;pars.xi],'%.5f'),...
                                                    num2str([sd;0],'%.5f'),...
                                                    strvcat(num2str(ObjAllMoms,'%.3f'),repmat(' ',[numel(Xi) 1])),...
                                                    strvcat(num2str(ObjTargMoms,'%.3f'),repmat(' ',[numel(Xi) 1])));
    Table6.Properties.VariableNames         = ["Point estimate" "sd" "Obj (all)" "Obj (targ)"];
    Table6.Properties.RowNames              = {Nu.targeted_pars{:},'xi'};
    disp(Table6)
    fprintf('\n')

    % Table 7
    Table7                                  = table(num2str(MomStats(:,1),'%.3f'),...
                                                    num2str(MomStats(:,2),'%.3f'),...
                                                    num2str(MomStats(:,3),'%.3f'),...
                                                    num2str(MomStats(:,4),'%.3f'),...
                                                    strvcat(num2str(OID,'%.3f'),repmat(' ',[Nmoments-1 1])),...
                                                    strvcat(num2str(dfr,'%.0f'),repmat(' ',[Nmoments-1 1])),...
                                                    strvcat(num2str(pval,'%.3f'),repmat(' ',[Nmoments-1 1])));
    Table7.Properties.VariableNames         = ["Emp. value" "Sim. value" "sd" "p-val" "OID" "dfr" "OID p-val"];
    Table7.Properties.RowNames              = Nu.all_moments;
    disp(Table7)
    fprintf('\n')

end

end